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Abstract. The Green Nagdhi equations are frequently used as a model of the 
wave-like behaviour of the free surface of a fluid, or the interface between two 
homogeneous fluids of differing densities. Here we show that their multilayer 
extension arises naturally from a framework based on the Euler Poincare theory 
under an ansatz of columnar motion. The framework also extends to the travelling 
wave solutions of the equations. We present numerical solutions of the travelling 
wave problem in a number of flow regimes. We find that the free surface and 
multilayer waves can exhibit intriguing differences compared to the results of 
single layer or rigid lid models. 
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1. Introduction 

Internal gravity waves have been observed propagating in many different locations 
in the world's oceans, both through direct measurement of the change in density 
stratification [Ij and through the observed change in surface conditions, from satellite 
observations [2jand from recent observations from the Space Shuttle in the region of 
Dongsha in the South China Sea [3]. These waves play a key role in the transport of 
energy in the ocean and, through wave breaking, help to control mixing [I]. The 
observed waves are not classical one dimensional phenomena. With wave crests 
extending up to 200 km perpendicular to the direction of motion, the wave properties 
vary with depth, exhibit curvature and interact with the underlying bathymetry and 
with each other. These interactions give rise to a number of phenomena including 
refraction, diffraction and wave front reconnection. 

This paper examines one strongly nonlinear, multilayer, two dimensional equation 
set for the behaviour of such waves, which is derivable from an Euler Poincare 
variational principle and shows that the equations admit travelling wave solutions 
which exhibit many interesting phenomena. We begin by discussing the derivation 
of the multilayer Green Nagdhi equations from a variational principle under the 
constraint of columnar motion. We then show the dynamics contain a strong, 
fast barotropic mode and relate this to the rigid lid models of other investigators. 
Finally we present numerical solutions to the travelling wave problem for the system, 
illustrating the important behaviour of the free surface. 

2. The multilayer Green Nagdhi equations 



We will begin by deriving a multilayer extension of the shallow water wave 
equations commonly attributed to Green and Nagdhi [5], using a variational principle 
under the ansatz of columnar motion. Consider a system of N homogeneous fluid 
layers which are acted upon by a constant gravitational acceleration, g and with the 
layers initially arranged in a stable stratification. Each layer thus possesses a constant 
density, pi for i = 1, . . . , N, with p\ < p2 < ■ ■ ■ < Pn- We have chosen by convention 
to number the layers downwards from the free surface (see Figure [TJ. Let (ui,Wi) 
denote the full three dimensional velocity fields within the layers and hi denote the 
depth of the layer interfaces with respect to the mean free-surface height at z = 0. 
We define b — — h^+i to be the fixed bottom topography to which the lowest layer is 



z{m) 





Figure 1. The multilayer fluid system 



An Euler Poincare framework for the multilayer Green Nagdhi equations 3 

assumed to remain attached. Summing over the layers, the total kinetic energy of the 
system is 

N r r r h, 



while the total gravitational potential energy, relative to a background state with 
vanishing density, is 

p p N phi p N 

/ Vdx := / J> / zdzdx = / E " ^+i) dx ' 

f=l J h i+ i J i= i Z 

We also assume two auxiliary equations within each layer, namely a three dimensional 
incompressibility condition, 

dwi , „ 

V • u, + — * = 0, 1 

az 

and a transport equation for the layer thicknesses, A = hi — hi+i, based upon the 
conservation of fluid within each layer, 
dD 

+ V ■ D iUi = 0. (2) 

We now make the ansatz of columnar motion with respect to the horizontal 
components of velocity in each of the layers, so that 

^ = fori = l,...,iV. (3) 

Together ([lj and ([3} enforce a linear dependence on z for the vertical component of 
velocity, Wi, so that using the bottom boundary condition, 

un ■ Vb + wn = 0, 

and integrating upwards we can find the vertical velocity within each layer in terms 
of the horizontal velocity divergences within that layer and those below it, 

N 

Wi = — zV • Ui — V • buj — V • Dj(uj — u,). (4) 

This gives a velocity jump across the interfaces of 

Wi(h i+ i) - w i+1 (h i+1 ) = (u i+1 - ■ V/ii+i, 

which is necessary to avoid the layers separating. Substituting ([3} into the definition 
of kinetic energy and integrating in the vertical gives 

/ Kdx= fj2 f (a|u 4 | 2 + A[^(V ■ u. ( ) 2 + AV ■ Ul W % + Wf]\ rfx, 

i— 1 ^ ^ 

where we have defined quantities 

(N \ N 

-b+ Dj + V-Dj-Uj, 
i=i+l / j=i+l 

which are equal to the vertical velocity at the bottom of each layer. 
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We now use Hamilton's principle in the form 
SS :=S / idxdt = 0, 



applied for the reduced Lagrangian defined by 

£ := K - V, 

coupled to The set of physically admissible horizontal velocities in each layer 
(a collection of vector fields over the horizontal domain, tangential on the boundary) 
form a Lie algebra under what is termed the ideal fluid bracket operator for given for 
functions F, G, of a vector field v^, by 

j- /--~f~i / \ [ fSG ^SF SF „SG\ , 

for variations defined by 

lim - [F(v + eSv) - F(v)\ = [ 5v ■ ^dx, 
e^o e J <5v 

while the layer thicknesses form a set of field densities which are Lie-transported 
(i.e. advected) by the Eulerian velocity flow. This places the system within the 
formalism of Euler-Poincare theory on semidirect products, which extends the results 
of Hamiltonian mechanics into this more generalized algebraic structure. 

Similar formulations have previously been derived for several alternative models 
for fluid flow, starting with the Euler equations, as originally considered by Poincare 
and including in particular the Camassa-Holm [6] and KdV [7] shallow water 
equations, models which, like the single layer Green Nagdhi equation, has been used 
successfully to model the behaviour of nonlinearly dispersive water waves. The Euler- 
Poincare theory has connections with many other topics in geometric mechanics, 
notably Lagrangian reduction [8], which has been applied successfully in several areas 
of compressible and incompressible fluid flow, and has implications for numerical 
methods, through a more rigorous understanding of the flow of conserved quantities 
in the system as geodesic motion on a suitable manifold. 

Calculating the requisite adjoint actions and their duals to generate the Euler- 
Poincare equations [9] in this framework (analogous to the Euler-Lagrange equations 
for a finite dimensional Hamiltonian system) gives the multilayer Green Nagdhi 
(MGN) equations of motion, 

dxrii _ T 5£ , , 

—^j- + Ui ■ Vm.; + riii • Vu t + niiV • = AV^, (6) 

in terms of the layer momenta, m.; := J^-, dual to the layer velocities. The MGN 
equation is forced by a pressure like term, -M^, containing the effects of gravity and 
nonhydrostatic terms, which appears out of viewing the relation between kinetic and 
potential energy in the system as a semi-direct product over the combined space of 
layer velocity and thickness fields. Explicitly, these terms are given by 

fD 3 WA \ 

Am - V f ■ Ui + A f -f^ ■ ^ + Wij X7h l+1 

3=1 V 
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i-l 



SDi 



i-i 

0=1 



Pi9h\ - ^2(pj - pi)Dj 
• u 



u* • V 



.7=1 



V 



D , 



V • u, + (uj - Ui ) 



Although complex, many of the terms in the equation are in balance. This is 
perhaps most clearly seen by rewriting (|6jl in the form 



d ( m t 
dt \ A 



m ■ v 



m, 

D, 



-D, 



Passing the material derivative part way through the definition of the layer momenta 
gives 

i-l 



d t Uj 
dt 



= -gV I hi + 



Pj ~ Pi 
Pi 



D 



j 



1 

Di 



-VDi 



dt 



di ( D, V ■ ii, 

dt 



111 
2 



Wi 



3 > 



identical to the equation set given for the same fluid system in the two layer case 
by Choi and Camassa [10J (henceforth the CC equation) by an asymptotic expansion 
method and by Liska and Wendroff \X1\ from directly substituting the definition of Wi 
in the Euler equations. It also agrees with the one dimensional multilayer equation 
of Choi [12]. The quantity labelled Ph is a hydrostatic pressure, representing forcing 
from both deviations to the free surface height and from variations in the thicknesses 
of the layers above. Calculating the full commutator [di/dt, Cijiij] for the symmetric 
second order elliptic operator defined by ir^ = C(D, b)ijUj, the equations may also be 



written entirely in terms of gradients of the horizontal velocity and layer thicknesses, 

1 Pj ~ P % | - 



C(D,b) 



diUi 
dt 




V D 



Ri 

T 



Si 
2 



sA Vh i+ i + vY,- D . 



R i 



where 



Ri := Di [(V • u,) 2 + tr(Vuf • Vu)] , 

a := - 



(7) 
(8) 



N 



• VV/i 1+1 ) • u, + V • (V • (DjUjUj)) 

i=i+i 

In this form we see the operator C acting as a smoother on the hydrostatic pressure, 
while there is an induced remainder term which grows with the non-linearity and 
further couples the layers. 
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A final rearrangement of (|6j gives an equation similar to the vorticity form of the 
Euler equations, 

51 



d ( mi \ _/xij • m; 
' - I + V 



V x 



m. 



XUi = V 5 Di 



dt\D l J \ D 
Taking curls shows that the equations materially conserve a quantity. 

1 /m, 

q, := — V x — — 
q Di \D t 

identified with potential vorticity in each layer. These conservation laws may also 
be derived from the existence of individual Kelvin circulation theorems in each layer, 
namely 

d f m, 

• dx = 0, 



dt Jc(u t ) Di 

where the integral is over the closed loop, c, assumed to move at the layer velocity, 
Uj. In turn it can be shown that the existance of these circulation theorems and their 
associated conservation laws follows in turn from the invariance of the EP formulation 
to fluid parcel relabelling in the configuration space, provided that it preserves Eulerian 
quantities. For fuller details of the Kelvin-Noether Theorem for EP equations in the 
context of semidirect product Lie algebras see [9]. The integral over the entire layer 
volume of any function of the qi represents a conserved quantity of the MGN equations. 
These are the Casmirs of the Lie-Poisson Hamiltonian operator, and hence of the Lie- 
Poisson bracket, ([5}. That is they are functionals, C, which satisfy {C,H} = for 
any Hamiltonian, H. 



3. The barotropic and baroclinic modes 



Linearizing the 2-layer MGN equations around a state of rest, with standing fluid 
heights di, di gives the system 



d_ 

dt 



1 , 1 



-gV(Di + D 2 ), 



d_ 



u 2 - -a^VV ' u 2 



P2 

dD 1 

^dt 
dD 2 
dt 



—d x { ^diVV ■ ui + d 2 VV ■ u 2 



c?iV • ui = 0, 



G? 2 V • u 2 = 0. 



-gV \ D 2 + ~Di ) , 



Taking divergences of the momentum evolution equations then give the linearized wave 
equations 

^ (Dx ^ADx ~ = gdxA(Dx + D 2 ), 
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Figure 2. The linear dispersion relation for the MGN equations. In this case 
di = 100m, pi = 1023kg m" 3 , d 2 = 200m, p 2 = 1026kg m" 3 . Note the wide 
separation of scales between the two modes. 



dt 2 



D 2 - 



rl 2 



/Mi .[di 



P-2 



-A -±D 1 + d a D l 



= gd 2 A D 2 



Pi 

P-2 



with A := d 2 /dx 2 + d 2 /dy 2 , the standard two dimensional Laplacian operator. 
Specializing to one-dimensional plane wave solutions, Di = Cie~xp[i(ujt — kx)], we 
obtain a linear dispersion relation 



1 + ^djk 2 



1 



dik 2 



P 1 J J 7„2 , 1j2j2»4 



Pi 



., d^ 2 k 2 + -d{d z 2 k 

3/92 9 \2p 2 



dfd^ LU 



1 



1 



d\d%k 



Pi 



2p 2 



di(di-d 2 yk' ! w 



-gk 2 ^(di + da) 
+gg'd 1 d 2 k 4 = 0, 

with </ := g(/9 2 — Pi) / p 2 a reduced gravity term. Plotting solutions shows both a fast 
barotropic mode, u>f, and a slow baroclinic one, u> s , as shown in Figure [2] Under the 
limits g'd\d 2 <C g(di + d 2 ) 2 , dik <C 1 we re-obtain the shallow water modes, 



I — y p- / d\d 2 

= V5(«i + d 2 jk, uj s = Jg'- 



-k. 



Looking for one dimensional solutions to the linearized system with u 2 = Xui and 
D 2 = (n — l)Di we obtain equations 



10 










3 2 






du i 




dx 



-did 2 



5^ 



»i 



+ X—did 2 

P2 



lei 

2p 2 



di 



-gp- 

dx 2 



8Di 

dx 



ui = {-gp + g') 



dDi 

dx 
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The final two equations give an obvious relation, 

(M-1)=^A, 
di 



o2 

while the first two equations can be rearranged to give in terms of quantities other 
than u, 



A 



(d? - 4) + A 



A _ pi 

2 P2 

2p 2 



dj 



d\d 2 

d d 2 ui 
dt dx 2 



= {gn{\-l)+g') 



8D X 

dx 



On re-substituting this gives a cubic equation for A. For g' small one approximate 
solution is /i = 1 + d 2 /d\,\ = 1, another consistent solution has /i 0, A » —di/d 2 . 
We identify the first solution with the barotropic mode and the second with the 
baroclinic. This shows the linear dynamics of the barotropic and baroclinic modes 
are governed by 



du bt 
dt ~~ 

and 

dLu bc 
dt 

respectively, where 



dD, 



-fl- 



it 



dx 

,dD b c 
' dx 



dD b 
dt 



= (di + d 2 ) 



du 



hi 



dx 



Lu = d\jd 2 



did 2 /3 



dD bc 
dt 

, Pi 



du bc 
dx 



—did 2 - ^—dl 

Pi 2/9 2 



d?_ 

dx 2 



u, 



the linearization of the MGN operator. This shows that to leading order the barotropic 
mode dynamics are precisely those of the one layer shallow water equations. Deviations 
in the free surface height are directly balanced by transient changes in the barotropic 
velocity. The baroclinic mode meanwhile shows the smoothing of the reduced gravity 
pressure, as represented in 

3.1. The rigid lid assumption 

As we have shown the free surface MGN equations exhibit a strong barotropic mode, 
which acts on a much faster scale than the baroclinic modes, which are the primary 
interest for studies of internal waves. This has important implications for numerical 
calculation, since for time stepping methods which are conditionally stable it will be 
the fast barotropic mode which sets the value of the constraint. For a kilometer deep 
ocean basin at hundred meter resolution Ax/ At ^ c = y/gH requires time steps below 
one second, much too small for most practical purposes. One method frequently used 
to modify the dynamics and remove this difficulty is to impose a rigid top boundary 
at z = 0. Under the EP framework this can be easily enforced using an additional 
constraint term in the reduced Lagrangian, 

' N 

^Dj-b dx.. 
,i=l 
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with <j)[x, y] a Lagrange multiplier be determined by the constraint that the term in 
brackets vanish everywhere . Under application of the EP methodology an additional 
pressure term equal to the multiplier cf>, thus constant in all layers, appears in the 
term due to the thickness densities, 5£/5Di. The value of <j> can be obtained through 
applying the rigid lid condition to the derived equations and solving the resultant 
elliptic system. The net effect of these additional terms is to modify the dispersion 
relation for the system, reducing its dimension so that the original fast mode is 
removed. The modified dispersion relation is found to be 



1 



1 



di + d 2 



- [d\d 2 


Pi _ 


1" 




.Pi 


H 3. 



d\d\ p\d\ 
2pT 



G 



9 



d\d 2 
di + d 2 



0, 



which possesses only the baroclinic roots, as stated. This means that the timestep for 
numerical methods is then controlled by the large baroclinic mode representing most 
of the energy of internal solitary waves. 

We observe that under the rigid lid assumption the no-normal flow condition 
requires that the vertical velocity vanishes on the top surface, so that 

w(z = 0) := £>iV • Ui + Wi = 0. 

If this substitution is made in the reduced Lagrangian before the variations are taken 
then the resulting equations in the two layer, one dimensional, case agree exactly with 
the equations from the multiplier method, since in one dimension conservation of the 
MGN PV is automatic. Both methods produce a system identical to the CC equation 
for flow under a rigid lid with varying topography. 



4. Travelling wave solutions of the MGN equations 

The one-dimensional GN equations are widely known to posses a travelling wave 
solution [5] of the form 



D = d 1 



— ; - 1 ) sech 2 
gd 



\/3(c 2 -gd) 



2cd 



(x — ct) 



1 



d 
D 



(9) 



where the quantity c is the group speed (and, since these are shallow water, waves 
phase speed) of a chosen wave. These waves exist for all c such that a condition for 
supercritical flow, 



gd 



> 1, 



(10) 



is satisfied. This also defines a Froude number, c/y/gd, for the system. The sech form 
is also found for the KdV equation, although the precise definition of velocity differs. 
The equation giving the GN travelling wave velocity, (JoJ) , follows directly from the 
transport equation for layer thickness Examining the Lagrangian (now calculated 
against a flat background state) in the case N = 1, 



/ 



I = / dxdy- 



D\u[ 



LP 
6 



du 
dx 



g{D-df 
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we see the travelling wave solution is a solution of the Hamiltonian formulation 
obtained by a direct substitution of the definition of u for a travelling wave, ([9]). 
The resulting Lagrangian is 

which under the standard Legende transform, p — gives a canonical Hamiltonian 
in terms of q — D of 



<ld2 -^-ch(i--\ - .7(7 -</r 



6<7 V 9, 

The same equation can also be derived indirectly [13] by noting that the one- 
dimensional GN equation possesses a conserved layer momentum, J m/Ddx, and that 
travelling waves are stationary functions of the quantity 

/m 
H — c—dx. 

That both methods give the same final result follows from the invariance of the original 
Lagrangian I to Galilean boosts, (x, t) — > (x—ct, t), so that changing to a frame moving 
at the constant wave speed requires only the redefinition of the rest kinetic energy. 

The Hamiltonian structure extends easily to the multilayer case, where each layer 
possesses an equation identical to ([9} . For general N the MGN travelling wave solution 
is a homoclinic orbit of the Hamiltonian system given by 



H(p,q) = -p i A- 1 p 



1 N 



N 



N 



,/ 1 - d j) 2 - J2 fa ^ d tf 

j=i j=i+l 



C 2 q 3 



1 _ it 

1] 



around the equilibrium point (p,q) = (0,d), where the p, q and d are the TV- 
dimensional vectors containing the pi, qi and di and A is an N x N matrix with 
elements defined by 



A, 



Pi d 

min{ij'} 



Ei— 1 
k=X 



k=l 



,2 

Pi~<h 
D k 

2D k 



* = J 

i 3 



Since A depends on q the system is not separable and the dynamics can be extremely 
complex. Figure [3] shows numerical solutions for the travelling wave problem, 
calculated by shooting along the unstable manifold from the equilibrium (p, q) = (0, d) 
at x = — oo. The integration problem was found to be extremely stiff and the use of a 
symplectic integration method, the generalized leapfrog method, was implemented to 
maintain stability. For comparison, we also plot the relevant numerical solutions to 
the CC equations [10] and the algebraic solution to the two- layer KdV equation |14] . 

The condition for supercritical flow and thus existence ( |Io| generalises to a 
condition that the stable manifold of the system at x = — oo must be of dimension 
greater than zero. This is equivalent to the condition that det(V — A 2 A) = has real 
solutions, where V is the matrix 



Vij = - 



d 2 n 



dqidqj 



(p,q) = (0,d) 



Pi {T z ~ 9 
9Pmin{i,j} 



i = j 

i + 3 
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Figure 3. Numerical shooting solutions for the MGN travelling wave problem 
showing (a) a two layer wave of depression (b) a two layer wave of elevation. The 
CC and two layer KdV solutions for the fluid interface are shown for reference. 
In both cases the interface of the MGN wave virtually overlies the CC wave 
for the chosen wave-speed with the KdV wave noticeably smaller in magnitude. 
However the MGN equations also give a free surface deviation. The apparent 
dipole structure of the free surface in case (b) is not evident in the individual 
layer thicknesses which remain precisely out of phase. 
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This matrix comes from the form of the Jacobian matrix of the Hamiltonian system 
linearized around the far field values. For N — 2 the critical condition for internal 
and external travelling waves become respectively 

cl_ > d 1 + d 2 T ((di + d 2 ) 2 - 4(p 2 - p 1 )d 1 d 2 /p 2 ) 1/2 

9 2 

This is however only a necessary condition on the wave speed and there exist 
stratifications for which no internal travelling wave is possible for c supercritical. This 
can be illustrated by considering the form of the potential 



1 N 

V(q)4£ 



JV 

Pi 



9 



j=i j=i+l 



1o 



in the Hamiltonian in the regimes where there exist waves of elevation, waves of 
depression, and for cases with no travelling wave. Contour plots of the potential 
in these three conditions are shown in Figure [4] along with the marked trajectories 
where a travelling wave exists. The potential is seen to undergo bifurcations with the 
creation or destruction of homoclinic contours around the equilibrium at q = d. The 
direction of the contour defines whether the wave is of elevation or depression and 
its destruction with increasing wave-speed represents a limit on the speed (and thus 
amplitude) of travelling waves allowed by the system. This follows since for N=2 the 
matrix A is positive definite and thus p T A _1 p is a strictly positive quantity. Jo and 
Choi [15] investigate the two-layer system in the rigid lid case and find conditions on 
minimum and maximum wave speed similar to those presented here. The critical flow 
condition is 

c 2 (p2 - p\)d x d 2 



9 P2(d 1 +d 2 ) 



this is the first order term in the Taylor series expansion of (111 in the limit 
(di + d 2 ) 2 ^> A(p 2 — p\)d\d 2 f p 2 . This suggests that the rigid lid model is a good 
approximation when density differences are small, or the aspect ratio d\jd 2 differs 
greatly from unity, as may be expected. A similar set of calculations and analysis is 
in progress for the case N = 3, which appears to increase the dimension of possible 
behaviour. 



5. Summary 

We have introduced a set of equations derived from a variational principle under an 
ansatz of columnar motion. These have been shown to be identical to the multilayer 
Green Nagdhi equations derived independently by other researchers by other methods. 
These equations are proved to contain a fast barotropic mode which is virtually 
unmodified by the nonlinear part of the MGN operator. This means the equations 
require careful handling in numerical simulations. The travelling wave solutions are 
shown to also be derivable from a variational principle, and to show a more complex 
range of behaviour and waveforms than the single layer or rigid lid cases. 
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Figure 4. Contour plots of potentials which allow (a) a wave of elevation, (b) 
a wave of depression (c) no travelling wave solutions. In the first two cases the 
trajectory for the travelling wave solution of the particular regime is plotted as a 
dotted line, contained within the zero contour. The fish-like looped zero contours 
disappear with increasing wave speed through a bifurcation with a second (hi, /12) 
contour joining (00, ±00) to (—00, ±00) with the sign of /12 positive for waves of 
elevation and negative for waves of depression. 
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